import os
import numpy as np
from skimage.feature import hog
from skimage.feature import local_binary_pattern
from skimage.feature import graycomatrix, graycoprops
from skimage.transform import resize
from skimage.measure import moments_central, moments_normalized, moments_hu 
from skimage.color import rgb2gray
from skimage import io

import pandas as pd
import shutil
import matplotlib as mpl
import matplotlib.pyplot as plt
import matplotlib.gridspec as gridspec
import matplotlib.colors as mcolors

from plot_t1t2_helper import plot_T1T2_distribution, add_box, add_right_cax, plot_T1T2

# 指定缩放后的大小
width = 64
height = 64
dim = (width, height)

def check_dir(file_path):
    if os.path.exists(file_path) ==False:
        os.mkdir(file_path)
    # else:
    #     # print("输出文件夹非空")
    #     shutil.rmtree(file_path)
    #     os.mkdir(file_path)


def read_excel_file(excel_file_Path):
    excel_data = pd.read_excel(excel_file_Path, sheet_name="T1-T2 relaxation spectrum",engine="openpyxl")
    return excel_data


def extract_hog_features(image):  
    # 加载图像
    
    gray = image
    # HOG特征提取
    hog_features, hog_image = hog(gray, orientations=9, pixels_per_cell=(8, 8),
                                cells_per_block=(2, 2), visualize=True, block_norm='L2-Hys')

    # print(hog_features.shape)
    # fig.show()
    return hog_features, hog_image

  
def extract_lbp_features(image):  
    # 加载图像
    
    gray = image
    # lbp特征提取
    radius = 3
    n_points = 8 * radius
    lbp_image = local_binary_pattern(gray, n_points, radius, method='uniform')

    # print(lbp_image.shape)

    return lbp_image


def visualization_features(image,features, features_name, file_save_path):
    # 可视化lbp特征
    fig = plt.figure(figsize=(12, 6),tight_layout=True)
    # 以10为底的对数，只显示正数
    log_axi_array = [0.01,0.1,1,10,100,1000,10000,100000]

    gs = gridspec.GridSpec(1, 2)
   
    ax_1  = fig.add_subplot(gs[0, 0])
    ax_2  = fig.add_subplot(gs[0, 1])

    # cax_2 = add_right_cax(ax_2, pad=0.02, width=0.02)
    # cax_1 = add_right_cax(ax_1, pad=0.02, width=0.02)
    cmp = plt.get_cmap('jet')
    #
    # # https://zhajiman.github.io/post/matplotlib_colorbar/
    #
    norm = mcolors.LogNorm(vmin=1E0, vmax=1E3)
    # # norm = mpl.colors.BoundaryNorm([0,1,2,3], cmp.N)
    # # create an Axes on the right side of ax. The width of cax will be 5%
    # # of ax and the padding between cax and ax will be fixed at 0.05 inch.


    # fig.colorbar(mpl.cm.ScalarMappable(norm=norm, cmap=cmp), cax=cax)

    im_1 = ax_1.imshow(image, cmap=cmp)
    ax_1.invert_yaxis()  #y轴反向
    ax_1.set_title('Original Image Feature')
    # ax_1.semilogx(base=10,subs=log_axi_array,nonpositive='mask')
    # ax_1.semilogy(base=10,subs=log_axi_array,nonpositive='mask')
    fig.colorbar(mpl.cm.ScalarMappable(norm=norm, cmap=cmp), ax=ax_1, label='Origin image colorbar')

    im_2 = ax_2.imshow(features, cmap=cmp)
    ax_2.invert_yaxis()  #y轴反向
    ax_2.set_title(features_name)
    # ax_2.semilogx(base=10,subs=log_axi_array,nonpositive='mask')
    # ax_2.semilogy(base=10,subs=log_axi_array,nonpositive='mask')
  
    fig.colorbar(mpl.cm.ScalarMappable(norm=norm, cmap=cmp), ax=ax_2, label=features_name +' colorbar')

    fig.savefig(file_save_path,dpi=220)
    # fig.show()
    return fig


# 同时显示hog和lbp特征
def visualization_features_joint(image,features_1, features_name_1, features_2, features_name_2, file_save_path,fig_dpi=220):
    # 可视化lbp特征
    fig = plt.figure(figsize=(18, 6),tight_layout=True)
    # 以10为底的对数，只显示正数
    log_axi_array = [0.01,0.1,1,10,100,1000,10000,100000]

    gs = gridspec.GridSpec(1, 3)
    
    ax_1  = fig.add_subplot(gs[0, 0])
    ax_2  = fig.add_subplot(gs[0, 1])
    ax_3  = fig.add_subplot(gs[0, 2])

    # cax_2 = add_right_cax(ax_2, pad=0.02, width=0.02)
    # cax_1 = add_right_cax(ax_1, pad=0.02, width=0.02)
    cmp = plt.get_cmap('jet')
    #
    # # https://zhajiman.github.io/post/matplotlib_colorbar/
    #
    norm = mcolors.LogNorm(vmin=1E0, vmax=1E3)
    # # norm = mpl.colors.BoundaryNorm([0,1,2,3], cmp.N)
    # # create an Axes on the right side of ax. The width of cax will be 5%
    # # of ax and the padding between cax and ax will be fixed at 0.05 inch.

    # fig.colorbar(mpl.cm.ScalarMappable(norm=norm, cmap=cmp), cax=cax)

    im_1 = ax_1.imshow(image, cmap=cmp)
    ax_1.invert_yaxis()  #y轴反向
    ax_1.set_title('Original Image Feature')
    # ax_1.semilogx(base=10,subs=log_axi_array,nonpositive='mask')
    # ax_1.semilogy(base=10,subs=log_axi_array,nonpositive='mask')
    fig.colorbar(mpl.cm.ScalarMappable(norm=norm, cmap=cmp), ax=ax_1, label='Origin image colorbar')

    im_2 = ax_2.imshow(features_1, cmap=cmp)
    ax_2.invert_yaxis()  #y轴反向
    ax_2.set_title(features_name_1)
    # ax_2.semilogx(base=10,subs=log_axi_array,nonpositive='mask')
    # ax_2.semilogy(base=10,subs=log_axi_array,nonpositive='mask')
    fig.colorbar(mpl.cm.ScalarMappable(norm=norm, cmap=cmp), ax=ax_2, label=features_name_1 +' colorbar')
    
    im_3 = ax_3.imshow(features_2)
    ax_3.invert_yaxis()  #y轴反向
    ax_3.set_title(features_name_2)
  
    # fig.colorbar(mpl.cm.ScalarMappable(norm=norm, cmap=cmp), ax=ax_3, label=features_name_2 +' colorbar')
    fig.colorbar(im_3, ax=ax_3, label=features_name_2 +' colorbar')

    fig.savefig(file_save_path,dpi=fig_dpi)
    # fig.show()
    return fig


# 矩阵翻转180度
def flip180(arr):
    new_arr = arr.reshape(arr.size)
    new_arr = new_arr[::-1]
    new_arr = new_arr.reshape(arr.shape)
    return new_arr

def flip90_left(arr):
    new_arr = np.transpose(arr)
    new_arr = new_arr[::-1]
    return new_arr

 

if __name__ == '__main__':
    
    headers = ["T1 time(ms)","T1 spectrum","T2 time(ms)","T2 spectrum"]
    orgin_data_path = os.path.join("塘油参1井-T1T2","T1T2-RELAXATION")
    # file_path_list: 存放原始项目数据路径下所有文件路径的结果列表
    file_path_list = []
    all_result_csvfile_name = ""  #结果文件
    # 项目数据第一步处理保存的位置
    features_fig_save_path = "features_figure"
    features_table_save_path = "features_table"
    orgin_fig_save_path = "T1T2Figure"
    check_dir(orgin_fig_save_path)
    check_dir(features_fig_save_path)
    check_dir(features_table_save_path)
    # 递归文件夹，找出excel的核磁数据文件
    count = 0
    nmr_dim = 64
    for parent, dirnames, filenames in os.walk(orgin_data_path,  followlinks=True):
        # 检测目标文件夹下的文件
        for filename in filenames:
            print("当前处理的文件是:",filename)
            # 读取核磁文件中对应的值
            nmr_data = read_excel_file(os.path.join(orgin_data_path,filename))
            T1_Times = nmr_data[headers[0]]
            T1_spectrum = nmr_data[headers[1]]
            T2_Times = nmr_data[headers[2]]
            T2_spectrum = nmr_data[headers[3]]
            nmr_dim = len(T2_spectrum)
            T1T2_data = nmr_data.iloc[:,len(headers):].values
            # print(T1T2_data)
            T1T2_fig, ax, ax_histx, ax_histy = plot_T1T2_distribution(T1_Times, T2_Times, T1T2_data)
            fig_name = filename.split(".xlsx")[0]
            T1T2_fig.savefig(os.path.join(orgin_fig_save_path,fig_name + "distribution.png"),dpi=220)
            T1T2_fig_na, ax = plot_T1T2(T1_Times, T2_Times, T1T2_data)

            
            fig_path_2 = os.path.join(orgin_fig_save_path,fig_name + "distribution_na.jpg")
            T1T2_fig_na.savefig(fig_path_2, bbox_inches='tight', pad_inches=0)
             
            # bbox_inches='tight', pad_inches=0 去白边
            # "./{}_resized.jpg".format(image_name.split(".jpg")[0]) 保留原始图像的名字并在后面加上_resized 

            # 原文链接：https://blog.csdn.net/qq_33377927/article/details/126512634
            
            # 方式一：直接使用原来的图
            T1T2_data_array = np.uint8(T1T2_data)
            
            # 方式二：读取matplotlib保存的图
            # T1T2_data_src = io.imread(fig_path_2)
            # T1T2_data_array_use = rgb2gray(T1T2_data_src)
            # T1T2_data_array = np.uint8(T1T2_data_array_use)
            # # T1T2_data_array_180 = flip180(T1T2_data_array)
            # T1T2_data_array_180 = np.flip(T1T2_data_array,axis=0)

            resized_img = resize(T1T2_data_array, dim, anti_aliasing=True)
            resized_img = T1T2_data_array
            print("resized_img.shape:",resized_img.shape)
            
            # 提取纹理特征
            glcm = graycomatrix(resized_img, [1], [0, np.pi/4, np.pi/2, 3*np.pi/4], levels=256, symmetric=True, normed=True)
            contrast = graycoprops(glcm, 'contrast').mean()
            energy = graycoprops(glcm, 'energy').mean()
            homogeneity = graycoprops(glcm, 'homogeneity').mean()
            dissimilarity = graycoprops(glcm, 'dissimilarity').mean()
            asm = graycoprops(glcm, 'ASM').mean()
            correlation = graycoprops(glcm, 'correlation').mean()

            # 提取结构特征
            mu = moments_central(T1T2_data_array)
            nu = moments_normalized(mu)
            hu_moments = moments_hu(nu)

            # 输出纹理+结构特征向量
            # feature_vector = np.concatenate((hist, [area, perimeter, aspect_ratio, contrast, energy, homogeneity,dissimilarity,asm], hu_moments))
            texture_structure_feature_vector = np.concatenate(([contrast, energy, homogeneity,dissimilarity,asm,correlation], hu_moments))
            # print(texture_structure_feature_vector)
            # print(texture_structure_feature_vector.shape)
            
            # 提取HOG特征
            hog_features, hog_image = extract_hog_features(resized_img)
            # visualization_features(hog_image, os.path.join(features_fig_save_path,fig_name + "_hog_image_" + ".png"))
                  
            # 提取LBP特征
            lbp_image = extract_lbp_features(resized_img)
            # visualization_features(lbp_image, os.path.join(features_fig_save_path,fig_name + "_lbp_image_" + ".png"))
            feature_figure_file_name = os.path.join(features_fig_save_path,fig_name + "_feature_" + ".png")
            
            fig = visualization_features_joint(resized_img,hog_image, "hog", lbp_image, "lbp", feature_figure_file_name,fig_dpi=300)
            label_text = "contrast, energy, homogeneity,dissimilarity,asm,correlation,Hus: "
            fig.text(60,10,label_text + str(texture_structure_feature_vector))
            
            col_names = ["contrast", "energy", "homogeneity","dissimilarity","asm","correlation","Hu_0","Hu_1","Hu_2","Hu_3","Hu_4","Hu_5","Hu_6"]
            pd_texture_structure_features = pd.DataFrame(texture_structure_feature_vector.reshape(1,-1),columns=col_names)
            
            pd_hog_feature = pd.DataFrame(hog_features)
            pd_hog_image = pd.DataFrame(hog_image)
            
            lbp_features = np.asarray(lbp_image).flatten()
            pd_lbp_feature = pd.DataFrame(lbp_features)
            pd_lbp_image = pd.DataFrame(lbp_image)
            
            # 分割结果
            # grayscale = rgb2gray(np.asarray(im_target_rgb))
            # pd_output_pred = pd.DataFrame(np.asarray(output_pred.detach().numpy()))
            # pd_im_target_gray = pd.DataFrame(grayscale)
            # pd_im_target_r = pd.DataFrame(im_target_rgb[:,:,0])
            # pd_im_target_g = pd.DataFrame(im_target_rgb[:,:,1])
            # pd_im_target_b = pd.DataFrame(im_target_rgb[:,:,2])
            
            writer = pd.ExcelWriter(os.path.join(features_table_save_path,fig_name + "_features.xlsx"))
            pd_hog_feature.to_excel(writer,sheet_name="hog_feature")
            pd_hog_image.to_excel(writer,sheet_name="hog_image")
            pd_lbp_feature.to_excel(writer,sheet_name="lbp_feature")
            pd_lbp_image.to_excel(writer,sheet_name="lbp_image")
            pd_texture_structure_features.to_excel(writer,sheet_name="pd_texture_structure_features")
            
            # pd_output_pred.to_excel(writer,sheet_name="output_pred")
            # pd_im_target_gray.to_excel(writer,sheet_name="im_target_rgb")
            # pd_im_target_r.to_excel(writer,sheet_name="im_target_r")
            # pd_im_target_g.to_excel(writer,sheet_name="im_target_g")
            # pd_im_target_b.to_excel(writer,sheet_name="im_target_b")
            
            writer._save()

            # break
            








